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ABSTRACT 

Two decades ago, empirical evidence concerning the existence and frequency of planets around stars, 
other than our own, was absent. Since this time, the detection of extrasolar planets from Jupiter-sized 
to most recently Earth-sized worlds has blossomed and we are finally able to shed light on the plurality 
of Earth-like, habitable planets in the cosmos. Extrasolar moons may also be frequent habitable worlds 
but their detection or even systematic pursuit remains lacking in the current literature. Here, we 
present a description of the first systematic search for extrasolar moons as part of a new observational 
project called "The Hunt for Exomoons with Kepler" (HEK). The HEK project distills the entire 
list of known transiting planet candidates found by Kepler (2326 at the time of writing) down to the 
most promising candidates for hosting a moon. Selected targets are fitted using a multimodal nested 
sampling algorithm coupled with a planet-with-moon light curve modelling routine. By comparing 
the Bayesian evidence of a planet-only model to that of a planet-with-moon, the detection process is 
handled in a Bayesian framework. In the case of null detections, upper limits derived from posteriors 
marginalised over the entire prior volume will be provided to inform the frequency of large moons 
around viable planetary hosts, 77^ . After discussing our methodologies for target selection, modelling, 
fitting and vetting, we provide two example analyses. 

Subject headings: planetary systems — planets and satellites: general — techniques: photometric — 
methods: data analysis — occultations 



1. INTRODUCTION 

Extrasolar moons represent an outstanding challenge 
in modern observational astronomy. Their detection and 
study would yield a revolution in the understanding of 
planet / moon formation and evolution, but perhaps most 
provoca tively they could be f requent seats for life in the 
Galaxy (Williams e t al.lll997h . Their existence is widely 
speculated both from a Cope rnican perspective and sim- 
ulations of planet synthesis (jElser et al.l 120111 ). and yet 
their detection has remained elusive. 

In truth, very little effort has been spent on the search 
for moons relative to their planetary counterparts, pre- 
sumably due to the known difficulty such a feat rep- 
resents. Indeed, there has never been a systematic 
search for exomoons and just a handful of papers have 
attempted to dete ct their presence in several transit- 
ing pla net systems dBrown et al.|[200lt Kipping fc~ Bakos 
1201 lallbT ) ([Lewis! 120111 provide a extensive discussion re- 
garding non-transiting moons around pulsar planets). 
With K epler now boasting 2326 candidate transiting 
planets (|Batalha et all f2012). we can consider this can- 
didate list to be comparable to target list of stars used 
in radial velocity surveys, for example. Further, Kepler 
is specifically designed to detect Earth-sized transiting 
objects and thus has the capability to find large moons 
(Kip ping et al.l 120091 ) . In light of this, the time is ripe 
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for the first systematic program to search for the satel- 
lites of extrasolar planets and so in this paper we describe 
our new observational project: "The Hunt for Exomoons 
with Kepler" (HEK). 

The objective of the HEK project is to inspect the most 
viable planetary candidates for evidence of an exomoon, 
using publicly available Kepler photometry. In cases of a 
null-detection, and when conditions permit, the obtained 
constraints will be used to inform a new statistic in exo- 
planetary science: the frequency of large moons, 77^ . 

The structure of the paper is as follows. In Sj2j we will 
present an overview of the current literature regarding 
the existence of exomoons, which will inform our search. 
In a similar vein, $3] presents an overview of the cur- 
rent literature on the observational consequences of ex- 
omoons. In Sj4j the multi-pronged HEK target selection 
process will be described, which greatly aids in refining 
the search to only the most favourable objects. $5] dis- 
cusses the detection, modelling and fitting strategy of 
the HEK project, with particular focus on the use of 
Bayesian evidence and the application of a multimodal 
nested sampling algorithm, MultiNest, in combination 
with our exomoon-modelling code, LUNA. In <|6l we ex- 
plore the process of vetting exomoon candidates against 
false-positives, such as perturbing planets. $7] presents 
two examples of the MultiNest algorithm with LUNA on 
synthetic data, illustrating the use of Bayesian evidence 
for detections. Finally, we summarise the project's goals 
and methods in Sj8] A list of acronyms used throughout 
this work is available in the Appendix, Table [3] Simi- 
larly, a list of mathematical symbols and notation used 
in this work is found in Table |4j 

2. EXOMOONS 
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2.1. Definition 

An extrasolar moon, or exomoon, is a natural satellite 
of a extrasolar planet. Our project will focus on the 
detection of satellites which are gravitationally bound to 
be within the Hill sphere of their host planet, as opposed 
to quasi-satellites which can reside outside. 

For a planet with a single large moon, it is possible 
that a more appropriate description would be a binary- 
planet. The blurry line between a binary-planet and 
a true planet-moon pair can be drawn by the point at 
which the centre of mass of two bodies lies outside the 
radius of both bodies. This distinction is merely a tax- 
onomical issue though and in principle the HEK project 
can also detect binary planets, should they exist. 

2.2. Observational Effects of an Exomoon 

A detailed discussion of how we will search for exo- 
moons is presented in f}3] However, for the purposes of 
providing some context in the subsequent sections, we 
will briefly outline the observational consequences of an 
exomoon. The HEK project will search for moons around 
transiting extrasolar planets and thus primarily make use 
of photometry (specifically photometry obtained by Ke- 
pler). For an o verview o f alterna tive exomoon detection 
techniques, see Kipping] (|2011b[ ). For transiting planet 
systems, a companion moon can reveal itself through two 
categories of observational effects: 

■ Dynamical variations of the host planet. 

■ Eclipse features induced by the moon. 

Dynamical effects are measured as perturbations of the 
motion of the host planet away from a simple Kcplerian 
orbit. It is thought that the most observable dynami- 
cal effects will be transit timing and duration variations 
(TTV and TP V respe c tively^ ([Sartoretti fc Schneiderl 
[1991 iSzabo et al.l [200l |Kj£rinj||2009a||bJ). Pynamical 
effects primarily reveal information about the exomoon 
mass. 

Eclipse features are caused by the moon either occult- 
ing the stellar light directly or inducing so-called "mu- 
tual events" {Ragozzine & Holman 2010) by occulting 
the planet during a planet-star eclipse. In either case, 
such events primarily reveal information about the exo- 
moon radius. 

2.3. Feasibly Detectable Systems 

Kip ping et al.l (|2009l ) conducted a feasibility study of 
Kepler's sensitivity to a habitable-zone gas giant with 
a single moon, based upon dynamical effects only (i.e. 
TTV & TPV) . The study assumed moons were on copla- 
nar, circular orbits around their host planet. They found 
that Kepler should be sensitive to exomoons of masses 
> 0.2 M®, in the most favourable circumstances. Even in 
this optimistic case, this is an order-of-magnitude more 
massive than the most massive moon in our Solar Sys- 
tem, specifically Ganymede at 0.025 Mg. Consequently, 
it must be understood that HEK is searching for moons 
in a mass regime which do not exist within our own So- 
lar System. We dub these moons as "large moons" and 
place a definition of such a moon as being > 10 _1 M®. 
Again, this is an arbitrary distinction but is useful for 
the HEK project. 



We point out that HEK will b e more sensitive than 
the limits in iKipping et al.l ()2009l ) due to the fact that 
we will also search for eclipses caused by the moon. This 
additional information means that we must be at least 
as sensitive as using dynamics alone, but should become 
significantly more sensitive by including this extra infor- 
mation. An in-depth analysis of our sensitivity to ex- 
omoons is beyond the scope of this work. We consider 
such an analysis unnecessary in light of the fact that 
each system studied will have upper limits derived and 
so the detectability of moons will become apparent as 
the project develops. At this time, we consider a far 
more useful application of our efforts to be in actually 
conducting the search. 

2.4. Plausible Origin of Large Moons 

Given that no large moons (Ms > lO^ 1 M®) exist in 
the Solar System, how likely is it that such an object 
exists? The first important point to make is that two 
classes of moons exist: regular satellites and irregular 
satellites. In discussing the plausibility of a large moon 
around a planet, there are two issues to consider: i) origin 
ii) evolution. In other words, the moon must have some 
initial plausible formation/origin scenario and second it 
must survive long enough to be detected. We here discuss 
the former issue. 

Regular satellites are those which are believed to have 
formed in-situ around the host planet a s it accumulates 
;as an d rock-ice solids from solar orbit. iCanup fc Ward! 
2006) predict that the cumulative moon mass is con- 
strained via J2 M S,t < 2 x 10~ 4 Af P , where M P is the 
mass of the planet host. This limit is caused by a balance 
of two competing processes; the supply of inflowing ma- 
terial to the satellites, and satellite loss through orbital 
decay driven by the gas. Thus, for Jupiter-like masses, 
large moons should not form. Although brown-dwarfs 
could harbour Earth-mass moons, the mass ratio remains 
too small for TTV/TPV methods to feasibly infer their 
presence. Therefore, we argue i t is unlikely HEK wil l 
detect any regular satellites if the ICanup fc Wardl (2006) 
scaling law holds. 

Irregular satellites are those which are obtained from 
a non-local origin suc h as capture (e.g. Triton; 
lAgnor fc Hami lton 2006j) or impact scenarios (e.g. the 
Moon: iTavlorl ll992jT Such moons could reach large 
masses so long as they are dynamically stable. This 
means that Earth-mass moons (or even larger) are plau- 
sible and HEK would be sensitive to such objects. Unlike 
regular satellites, irregul ar moons may frequentl y reside 
in retrograde orbits too (|Porter fc Grundvll2011l ). 

In order for a planet to capture a large moon, an 
essentially terrestrial-mass object must be dynamically 
captured by a larger body. The probability of such an 
event occurring is a topic of active theoretical research 
and HEK will elevate it to a topic of active observa- 
tional research too. Several plausible capture mecha- 
nisms have been pr o posed . From a case study of Triton, 
lAgnor fc Hamilton (|2006| ) propose that this moon was 
originally a member of a binary which encountered Nep- 
tune during its migration through the proto-Kuiper belt. 
Upon encountering Neptune, a momentum-exchange re- 
action occurred ejecting one member and bounding the 
other as a satellite. One difficulty with this mechanism 
is that a binary terrestrial planet pair is required to en- 
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dow an Earth-mass moon around a gas giant and the 
abundance of such binaries is unclear. 

Another possible capture process is via atmospheric 
tunneling, where a terrestrial mass object encounters the 
atmosphere of the gas giant, inducing str ong drag effects 
leadi ng to large changes in momentum (IWilliams et al.1 
1201 The obvious extreme limit of this scenario is an 
impact between two rocky cores, where the smaller mass 
body is broken up and later reforms in a close-in orbit, 
as is proposed for the Moon's origin (|Tavlorl H992l Im- 
pacts do seem to b e a feasible way of forming moons with 
lElser et"aT1 ()2011l ) recently simulating the interaction of 
planetesimals and finding that Earth-Moon pairs should 
be common. 

For planets which do not migrate through a proto- 
Kuiper belt or under the assumption that such objects 
will never reach sufficient mass to qualify as large moons, 
an alternative source of terrestrial mass objects is re- 
quired. This object could be an inner terrestrial planet 
encountered during the gas giant's inward migration or 
even a large, unsta ble Trojan which lib rates too close to 
the planet. Indeed, lEberle et al.1 (|201Q[ ) have shown that 
a gas giant planet (in their case HD 23079b) can cap- 
ture an Earth-mass Trojan into a stable satellite orbit, 
occurring in 1 out of the 37 simulations they ran. 

Any capture process (as opposed to an impact) will 
tend to produce very loosely-bound initial orbits, leading 
us to question how many of these captured bodies may 
ultimately survive. iPorter fc Grundvl (|2011h have inves- 
tigated this issue and found that captured moons have 
encouraging survival rates, with a survival probability of 
order 50% for various planet-moon-star configurations. 

Fi nally, true binary pla n ets ar e in principle also plausi- 
ble. iPodsiadlowski et al.l (|2010l ) showed that one viable 
scattering history in the formation of a planetary sys- 
tem is the tidal capture of two planets forming a binary. 
Indeed, a Jupiter-Earth pair could be considered as an 
extreme binary, much like Pluto-Charon. 

2.5. Plausible Evolution of Large Moons 

With several plausible avenues for a large moon to 
become bound to a planet, the next requirement is 
that the moon can surv ive long enough to be observed. 
IPorter fc Grundvl (2011) show that captured moons tend 
to start in eccentric, inclined orbits and rapidly circu- 
larise and relax into a coplanar orbit. With a moon on a 
circular, coplanar orbit within the Hill sphere, the moon 
evolves by either spinning-in or out through tidal effectfQ. 
If the rotational period of the planet is shorter than the 
moon's orbital period, the tides cause the moon to spin- 
out over ti me. The opposite proces s occurs if this ratio 
is reversed (jBarnes fc O'Brienl[2002l ). 

Regardless as to whether the moon moves inwards or 
outwards, the spatial boundary conditions must be given 
by the moon being in contact with the planet (minimum 
spatial separation) to being outside the sphere of grav- 
itational influence of the planet (the maximum stable 
separation). The time to move between these two lim- 
its is insensitive to the direc tion of the moon's migra- 
tion (jBarnes fc O'Brienll2002|) . A fortuitous moon could 
in principle reverse the direction of its migration just 

7 Tides can also induce significant heating on the satellite 
IjCassidv et al.ll2009h 



before hitting one of these spatial limits (due to brak- 
ing of the planet's rotation for example) and thus essen- 
tially double its maximum allowed lifetime. However, we 
do not consider such a scenario to be likely. Since the 
tidal dissipation depends upon the mass of the moon, 
one can write the maximum allowed exomoon mass as 
a function of the moon' s lifetime (expression taken from 
iBarnes fc O'Brienl l2002). This yields 
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where Af* is the stellar mass, Mp is the planetary mass, 
Qp is the tidal quality factor of the planet, lt2 P is the Love 
number of the planet, G is the gravitational constant, Rp 
is the planetary radius, as* is the orbital semi-major 
axis of the planet-moon barycentre around the host star 
and T is the lifetime of the moon. Equation [1] suggests 
that Earth-like (i.e. habitable-zone) moons are plausible 
around Jupiter s for billions of years aro und stars of mass 
M, > 0.4 M s (Ba rnes fc OTirienl 12001) . 

In Equation [TJ S) max represents the maximum stable 
semi-major axis for the moon, in units of the Hill ra- 
dius. One may naively expect this should be equal to 
unity but in reality three-body perturbations tend to dis- 
rupt a moon before it reaches t he Hill radius. Through 
detailed numerical integrations, iDomingos et al.l (2006) 
have shown that £> max = 0.4895 for prograde satellites 
and £> max = 0.9309 for retrograde satellites (assuming 
circular orbits). 

Finally a planet which migrates inwards will eventu- 
ally lose its moon(s) due to the shrinking Hill sphere. 
The Hill radius is given by 



Rh = «s* 



M P 
3M~, 
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(2) 



One can therefore see that the Hill radius decreases lin- 
early with the o rbital semi-major axis of the host planet. 
iNamounll (|2010l) showed that since planetary migration 
occurs much faster than moon migration, a moon ini- 
tially w ell-inside the Hi ll radius can quickly find itself 
outside. INamounl (|2010f ) estimate moons tend to be lost 
by this process for as* < 0.1 AU. 

The planetary parameters clearly have a significant im- 
pact on the survivability of a large moon. Since in general 
the planetary parameters may be reasonably estimated 
based upon the transit light curve, it is possible to create 
a list of the most favourable planetary candidates for ex- 
omoon inspection. More details on this target selection 
procedure are given in Sj4] 

2.6. Objectives of HEK 

The existence of large moons is hypothetically plausi- 
ble, but currently we have no empirical evidence to test 
this hypothesis. For this reason, the objectives of the 
HEK project will be as follows: 

1. The primary objective of HEK is to search for sig- 
natures of extrasolar moons in transiting systems. 
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2. The secondary objective of HEK will be to derive 
posterior distributions, marginalised over the en- 
tire prior volume, for a putative exomoon's mass 
and radius, which may be used to place upper lim- 
its on such terms (where conditions permit such a 
deduction) . 

3. The tertiary objective of HEK is to determine 77^ - 
the frequency of large moons bound to the Kepler 
planetary candidates which could feasibly host such 
an object (in an analogous manner to 77^ - the fre- 
quency of Earth-like planets). 

For the primary objective, the issue of upper limits 
or detection biases is irrelevant and in many ways this 
is the simplest task to execute on candidate systems. 
In contrast, the secondary objective requires more care 
due to the inter-paramet er correlations. For example, 
iKipping fc Bakosf (|2011b[) show how the excluded 3-cr 
limit on a putative moon around TrES-2b is strongly cor- 
related to the assumed semi-major axis and inclination 
of the moon. For this reason, any upper limits must re- 
late a posterior marginalised over the entire prior volume 
(more details on this are given in i)5.3l) . Additionally, in 
some cases it may not be possible to provide a mass or 
radius upper limit if, for example, a strong eclipse signal 
is found but ultimately deemed to be a false positive. 
Finally, the tertiary objective is challenging in light of 
the numerous detection biases which plague any survey 
such as this. For example, two problematic biases come 
from the fact we preferentially select systems with short- 
cadence data (see ^4.6[) and visual anomalies (see £14.31) 
in the light curve. 

In the case of a definitive signal, we require a detection 
significance thresho ld. Currently, 2326 planetary can- 
didates are known (Batalh a et al.ll2~012[ ). but this num- 
ber is constantly increasing. We may conservatively set 
our false- alarm-probability threshold to be 1 in 10,000 or 
3.89-er. This is rounded up to 4-cr for our nominal detec- 
tion threshold. Statistically, this implies that 1 in every 
15,787 claimed detections will be false, although in real- 
ity the false-positive-rate is a function of how carefully 
we vet our candidate signals rather than merely the num- 
ber of sigmas the signal is detected to. Acknowledging 
this, we will follo w the set of detection criteria defined in 
Kipping (2011b). Most relevant of these, the exomoon 
targets must be physically feasible solutions (criterion 
C3). Details on our vetting process, which will undoubt- 
edly evolve as the HEK project develops, are given in 

m 

In the case of a null-detection, the objective of HEK 
will be determine the exomoon mass and radius which 
can be excluded, which will aid in the determination of 
r/(j . As already mentioned, these upper limits will be 
based on a posterior marginalised over the entire prior 
volume. Using a 4-cr upper limit is usually not very 
informative since they allow for "hidden moons" , for 
example a moon which is behind the planet in every 
transit epoch (e.g. see §7.1j) . In such a case, the ra- 
dius limit is essentially unbound. We therefore choose 
to give the lower constraint of 90% confidence upper 
limits. Although less stringent, these limits are more 
useful in a population-sample of exomoon candidates. 
Further, all null-detections will have the posterior dis- 



tributions made publicly available on the project web- 
site (www.cfa.harvard.edu/HEK/) so that the commu- 
nity may investigate rjg using our results. 

2.7. Why HEK is Possible 

HEK is feasible for two reasons. Firstly, as already 
discussed, Kepler has yielded extraordinary success with 
2326 transiting cand idates down to Earth-sized objects 
(|Batalha et al.ll2~012D . Secondly, recent advancements in 
the theoretical development of exomoon-search methods 
make a search feasible with modern co mputational re- 
sources. Specifically, the LUNA algorithm (lKippingl2011a| ) 
offers a completely analytic and exact solution for the 
planet-moon light curve including dynamics, non-linear 
limb darkening and the modelling of mutual events. 
Methods based upon even partial numerical implemen- 
tation (such as pixelating the star) dramatically impinge 
our ability to run light curve fits, given the large number 
of parameters and prior volume which must be explored 
in any given fit. Thus, Kepler and LUNA are both critical 
to the feasibility of HEK. 

3. OBSERVATIONAL EFFECTS 

3.1. Dynamical Variations 

As discussed in ^2 -21 there exists two broad categories 
of observational consequences of an exomoon in the tran- 
sit light curve: dynamical variations and eclipse features. 
In this section, we will overview these techniques, each 
of which is a key tool to the HEK project. We begin our 
discussion with dynamical variations, of which there are 
several flavours. 

3.1.1. Transit timing variations (TTV) 

The first conceived effect was transit ti ming variations 
(TTV), bv iSartoretti fc Schneider] ()1999D . which is con- 
ceptually analogous to the astrometric technique of find- 
ing planets around stars. The motion of the planet 
around the planet-moon barycentre causes a planet to 
transit perio dically early and late. For co -aligned, cir- 
cular orbits, ISartoretti fc Schneiderl (|1999f ) showed that 
the maximum deviation in the time between two transits 
would be 

~ 3V37rMp /3 M* 1/3 

V M e / V years / V M P ) \ M* / 

(3) 

where Pb* is the orbital period of the planet-moon 
baycentre around the star (the "B" subscript will be 
used to refer to the planet-moon barycentre through- 
out). Encouragingly, planets on orbital periods of ~ 1 
year coul d genera t e very significant TTV amplitudes. 
However, IKipping] (|2009af) pointed out two critical hur- 
dles with the technique. Firstly, many other effects can 
cau se TTVs aside from moons, notably pertu rbing plan- 
ets (|Agol et alj2005HHolmari fc Murravl2005D . and there 
seemed to be no obvious way to discriminate between 
the physical source of the measured TTV. Secondly, the 
satellite's period around the barycentre {Psb) is less than 
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the p lanet's orbital p eriod for all bound orbits. Specifi- 
cally, Kipping ( 20093) have shown that 



Psb 



/©3 



(4) 



since J) < 1 for all bound exomoons, the period of the 
exomoon will always be less than 60% of the planet's 
period. This means the TTV waveform is undersam- 
plec0 since one can only measure a timing deviation once 
per transit. This undersampling means that one can- 
not reliably infer the period of the exomoon signal (only 
a set of harmonic solutions) and instead one can only 
reliably measure the RMS amplitude (i.e. the scatter) 
which scales as ~ cisbMs- Whilst knowing ass Ms is 
useful, clearly knowledge of each component would be 
more powe rful in und erstanding an exomoon. 

IKippinel (1201 lbl) gene ralized the original 
ISartoretti &: Schneider! (|1999D equation for the am- 
plitude of the TTV effect, to account for longitude of the 
ascending node, in clinations, eccen tricities and position 
of pericentres (see [Kippingj 1201 lbl for definitions of the 
various terms). The root-mean-square (RMS) amplitude 
of the full TTV effect is given by 



3 TTV 



clsbMsPb* (1 ~ ejgV 1 - 4. 
clb*M p (1 + e B * sinujB*) 
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, (5) 



where the $ term contai ns information on the eccen- 
tricity and can be found in IKippind (|2011bD . The asso- 
ciated waveform is given by 



TTV: 
where 
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A T tv = [1 + esscos^ss] _1 [ cosn7 ss cos(v S b + u S b) 

- smi S BsmzusBsm(v S B + ^sb)}- (7) 

where i>sb is the true anomaly of the satellite around 
the planet-moon barycentre. 

3.1.2. Velocity induced transit duration variations (TDV-V) 

TDV-V was c onceived of a d ecade after TTV by 
Kipping (2009a). Kipping (2009a) showed that the same 
motion responsible for changes in position causing TTV 
should also cause changes in velocity. Since the planet's 
velocity is inversely proportional to the duration of the 
transit, TDV-V must be another observational conse- 
quence of extrasolar moons. Whereas TTV scales as 
cisbMs, thus favouring the detection of moons at large 

— 1/2 

separation, TDV-V scales as a SB Ms and so favours 
the detection of close-in moons. TDV-V is conceptually 
analogous to the radial velocity method of finding planets 
(except really it is tangential velocity). 

Besides from the complementarity of their parame- 
ter space coverage, TTV and TDV-V could also yield 

8 This applies to not only TTV, but for all exomoon-induced 
timing effects i.e. TDV-V and TDV-TIP 



a unique solution for M$ and asB, if both effects were 
detected 0- This therefore solves the undersampling is- 
sue presented by only detecting one of the two effects, 
as discussed earlier for TTV. Further more, for coplanar, 
circular moons, the two effects exhibit a phase difference 
of 7r/2 radians since the dynamics is essentially projected 
simple harmonic motion. This phase difference offers a 
method to distinguish an exomoon TTV from, say, a per- 
turbing planet induced TTV. Therefore, detecting both 
TTV and TDV-V solves the issue of solution unique- 
ness (assuming the moon is coplanar and circular) and 
mass determination. Defining the Tb as the duration 
for the planet-moon barycentre to enter then exit the 
stellar disc, the TDV-V waveform may be described via 
(Kip pin3l20TTbl) to be 



TDV - V = T E 



clsbM s Pb* 
a,B*MpPsB 



v/T^ 



V 1 ~ e SB C 1 + e B* sinw B *) / 



A 



TDV-V, 



(8) 



where 



Atdv-v = cosnjsslesssinwss + sm{v SB +usb)} 

+ smi S B smw S B[esB coswss + cos(v S b +u> SB )} 

(9) 

The correponding RMS ampl itude may be found 
through integration over v$b (see [Kip ping 20 1 lbl for de- 
tails) and yields 
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where the $ term again cont ains information on the 
eccentricity and can be found in IKippind (|2011b[ ). 

3.2. Transit impact parameter induced transit duration 
variations (TDV-TIP) 

Later, [Kipping (2009b) showed that a second tran- 
sit duration variation (TDV) effect exists, dubbed 
TDV-TIP, standing for transit-impact-parameter in- 
duced TDV. The physical origin of this effect is that if the 
planet-moon orbital plane is not precisely normal to the 
sky (which is practically always true), then the planet's 
reflex motion must yield a component orthogonal to both 
the observer's line-of-sight and the tangent to the plane- 
tary motion. This motion leads to periodic variations in 
the apparent transit impact parameter. Since the tran- 
sit duration is a strong function o f the impact parame- 
ter (jSeager k. Mallen-Ornelasl2003l ). then duration varia- 
tions must follow. The TDV-TIP effect is typically much 

9 Psb may then may easily calculated through Kepler's Third 
Law 
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Fig. 1. — Cartoon of the TDV-TIP effect. Here, the moon is re- 
laxed into the same orbital plane as the planet's orbit and causes 
the planet to experience reflex motion illustrated by the two posi- 
tions of the planet, (1) and (2). This motion can be seen to cause 
a change in the apparent impact parameter, which causes a change 
in the transit duration. 



smaller than the TDV-V (velocity) effect but is strongly 
enhanced for near-grazing transits. Intrigueingly, TDV- 
TIP allows one to determine the sense of orbital motion 
of the moon (i.e. prograde or retrograde) and thus could 
be a key tool in understanding the origin of the satellite 
(although such cases require high signal-to- noise) . The 
phase of TDV-TIP is constructive to TDV-V for pro- 
grade orbits and destructi ve for retrograde . The TDV- 
TIP waveform is shown in Kipping (20lT[J) to be 



TDV - TIP 
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where the $ term again cont ains information on the 
eccentricity and can be found in lKippind (|2011b[ ). 

It should noted that due to the undersampling issue 
discussed in regard to TTV, TDV-V and TDV-TIP are 
simply observed as a global TDV effect. Thus we have 
two timing observables, TTV and TDV. 

3.3. Eclipses Features 

A second class of observational effect we can search for 
is the eclipse of the moon. This eclipse can be either 
in- front of the star or in-front/behind the planet during 
the planet-star transit. The former, which we refer to 
as "auxiliary transits" , is more likely to be detected for 



moons on wide orbits. The latter, sometimes dubbe d 
"mutual events" (|Ragozzine fc Holman|[20i0t |Pall2011ft . 
is geometrically more probable for moons on close-in or- 
bits. We show examples of each of these types in Fig- 
ured 

Critically, unlike TTV and TDV effects, eclipses are 
sensitive to the exomoon radius (and not the mass). 
This presents the opportunity of determining the density 
of the moon. Further, in an analogous manner to how 
transiting planets may yield the mean stellar density, p* , 
(jSeager fc Mallcn-Ornclas 2003), transiting moons allow 
one to measure t he mean planeta ry density directly, pp, 
as first shown in iKippingj (|2010cD . Armed with both p* 
and pp plus the ratio of the -Rp/i?*, one can determine 
Mp/M* (|Kipping|l2010d ). This powerful trick means that 
even if no TTV or TDV signal is detected, the eclipses of 
the moon alone allow one to measure the mass of the exo- 
planet (assuming M* is known). If radial velocities exist 
for the system, can be me asured directly through 
this technique (|Kipping| l2010cD . Additionally, the dy- 
namically determined Mp/M* can be compared to the 
RV solution using a-priori value of AT* (say from spec- 
troscopy) to ensure consistency. 

The ability to weigh planets using moons plays a key 
role in the HEK project. Candidates with unphysical 
properties can be quickly rejected as plausible exomoon 
signals. 

TIP, 4. TARGET SELECTION 

4.1. Candidate vs Planet 

At the timing of writing, 2326 Ke pler transiting cand i- 
date planets have been announced (jBatalha et al.l[2012f ). 
However, the p ublic candidate list is only 1235 at the 
time of writing (|Borucki et al.l [201 lh . HEK will employ 
the most up-to-date list during the development of the 
project. It should also be noted that included within the 
candidate list are several candidates which have already 
been confirmed as bona-fide planets to a high degree o f 
confidence e.g. the Kepler-9 system (|Torres et alJl20lTl ). 

Despite the progress made in confirming these ob- 
j ects through nov el techniques, such as BLENDER 
(|Torres et al.l [20041. the vast majority of the planetary 
candidates remain unconfirmed. 

One possible approach might therefore be to only in- 
spect the confirmed planets for exomoons. However, such 
a strategy is unnecessarily restrictive. The detection of 
an exomoon would in fact allow one to measure the mass 
and ra dius of the planet, as well as the moon ((Kipping 
l2010d ). By modelling these systems, one essentially ob- 
tains the same information which is normally provided 
by radial velocity (RV) or multi-planet transit timing 
variations (TTV) - the massed. Therefore, exomoon 
detection is a planetary confirmation tool, in addition 
to the RV and TTV techniques. Consequently, there is 
no need to limit ourselves to confirmed planets only and 
the HEK project will consider all planetary candidates 
as possible exomoon hosts. 

4.2. Target Selection (TS) Overview 

Equipped with the capacity of exomoons to confirm 
planetary candidates, analyzing all of the hundreds of 

10 Stictly speaking, one obtains the mass ratios 
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Fig. 2. — Left panel: transit light curves of a planet with a moon on a wide separat ion, demonstrati ng auxiliary transits. Right panel- 
transit light curves of a planet with a close-in moon, demonstrating mutual events. See \Kii)'Dind f 20 lid) for details of the parameters used 
in these simulations (Figures 7&6 respectively). 



planetary candidates for evidence of exomoons would be 
the most comprehensive way to proceed. However, prac- 
tical limitations, such as man-power and computational 
constraints, make such a task unrealistic. Therefore, we 
must conduct a target selection (TS) phase. Target se- 
lection works by taking the full list of pl anetary candi- 
dates (in our case the KOIs presented in iBorucki et al.l 
1201 It ) and identifying those targets which are of highest 
priority for more intensive investigation. We have three 
strategies for target selection, which have some overlap: 

1. Visual inspection a subset of planetary candidates 
for exomoon-like eclipse features (TSV). 

2. Automatic filtering of the planetary candidates, 
based upon system parameters (TSA). 

3. Targets of opportunity (TSO). 

In some cases, a planetary candidate may fall into 
more than one category. The identified candidates are 
then further prioritized by hand. This final stage, which 
we call target selection prioritization (TSP), is typically 
done by experience of what constitutes a feasible candi- 
date (see £|4.6I for more information). 

4.3. Visual Target Selection (TSV) Overview 

Visual Target Selection (TSV) is typically conducted 
by first selecting sub-sample of planetary candidates. 
This subset is typically selected from some simple pa- 
rameter constraints such as planet size and semi-major 
axis. The aim is to allow this subset to have only a weak 
theoretical prejudice for where to expect a moon. The 
subset sample size may range from dozens to hundreds 



of light curves. The light curves are then inspected for 
features resembling exomoon-like signals, notably: 

■ Mutual events, preferably with a flat-top (to dis- 
criminate from a starspot crossing) e.g. right-panel 
of Figure H 

■ Auxiliary transits (a second distinct transit feature 
offset from the primary) e.g. left-panel of Figure [2j 

■ Repeating anomalies in the light curve of any mor- 
phology. 

The visual inspection of Kepler photometry has been 
already been shown to be a useful technique fo r find- 
ing transiting pla net candidates (|Fischer et al.l 120121 : 
iLintott et aLl l2012) and here we extend such inspections 
to hunting for moons too. Whereas TSA relies heavily on 
the derived parameters of the system (plus an assumed 
mass-radius relation), TSV is far more ignorant of these 
values. The benefit of this is that we are able to catch in- 
teresting systems which would go missed by TSA due to 
potential errors in the assumed system parameters (for 
example from the Kepler Input Catalogue, KIC) or even 
our theoretical understanding of where moons may re- 
side. TSV systems are later scrutinized to see whether 
the signal is sufficiently high signal-to-noise, in particular 
relative to any time correlated noise in the photometry, 
but this forms part of our final prioritization stage, TSP 
(see £14.61 for more information). 

4.4. Automatic Target Selection (TSA) Overview 

TSA works by applying a set of filters which reduce the 
total planetary candidate list to to a more manageable 
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size. There are several considerations affecting the target 
selection: 

■ Availability - Systems for which there exists avail- 
able photometry. 

■ Reliability - Systems for which we have reasonable 
parameters. 

■ Capability - Systems which are capable of hosting 
large moons. 

■ Detectability - Systems which are capable of pre- 
senting a detectable moon signal. 

TSA leans heavily on the system parameters. Notably, 
the planetary mass must be essentially guessed based 
upon the radius. Details on this procedure are presented 
in §4.4.2. 

Stellar parameters of the planetary candidates have 
been estimated by the Kepler team as part of the Kepler 
Input Catalogue (KIC), a photometric survey of stars 
in the Kepler field-of-view designe d to identify bright 
dwarfs as targets for the mission. iBrown et all (|2011[ ) 
describes the estimation o f physical paramete r s, whe reby 
the synthetic spectra of iCastelli fc Kuruczl (|2004fl are 
forward-modelled with effective temperature, T c g, sur- 
face gravity, log(g), and metallicity, log(Z), as free pa- 
rameters to match the photometric measurements in the 
KIC. A relation between luminosity, effective tempera- 
ture and surface gravity, deriv ed from the s t ellar e vo- 
lutionary models computed by iGirardi et al.l ()200ClO . is 
used to estimate the stellar masses, which, when com- 
b ined with log ( g), est imates the radii. 

IBrown et al.1 (|2011[ ) state that the KIC effective tem- 
perature and radius estimations are reliable for Sun-like 
stars, but are "untrustworthy" for stars with T c g < 
3750 K. To accommodate this possible weakness, TSA 
is run on both the KIC-only catalogue a nd a catalogue 
of KIC plus that of lMuirhead et al.l (|20111 ) , who use near- 
infrared spectroscopy. These cases will be flagged appro- 
priately. By repeating for both catalogues, we make no 
decision about which catalogue is the correct one and en- 
sure we do not miss any interesting objects. Any other 
catalogues which appear will be treated in a similar man- 
ner. 

4.4.1. Availability 

At the time of writing, the only publicly available pho- 
tometry from Kepler comes from quarters 0, 1, 2 & 3 
(Q0, Ql, Q2 & Q3). In all cases, Ql (33.5d), Q2 (88. 7d) 
and Q3 (89.3d) photometry exists and naturally this will 
continue to expand as time goes on. Only in some cases 
is Q0 (9.3d) photometry available. This is because Q0 
was technically "calibration" data and originally not ex- 
pected to be useful scientifically. As a result, the number 
of sources observed was only 52,496 in Q0 but expanded 
to 156,097 for Ql. Whilst some sources were dropped 
from the target list, the majority were kept and thus 
roughly a third of all sources have Q0 photometry avail- 
able, in addition to the other quarters. 

In most instances, photometry data is currently avail- 
able in long-cadence (LC) mode only, whereas SC is 
preferable for our search for reasons described in ^4.61 
There are three possible exceptions to this rule: 



1. A planetary candidate was detected very early on 
and thus the data was able to be switched to short- 
cadence (SC) for subsequent observations. 

2. The star already had a known exoplanet and thus 
was observed in SC from the outset. 

3. The star was coincidentally pre-selected for astero- 
seismology and thus was observed in SC for one or 
more of the three available quarters. 

Case 3) is not very useful to us because the subset is 
rather small. Case 2) is also not useful because the known 
exoplanets, HAT-P-7b, HAT-ll-b and TrES-2b, are all 
short-period hot-Jupiters which are unlikely to be good 
exomoon hosts. Case 1) does not guarantee SC data be- 
cause only 512 targets can be observed in this mode at 
any one time and the number of pla netary candidates ex - 
ceeds this value significantly f2326; lBatalha et alj [2012). 

In conclusion, it is generally unlikely that a potentially 
interesting target (for exomoon hunting) will have any 
useful SC data available at the start of the HEK project 
and this represents a major limitation on our capabili- 
ties. Those that do are highly valuable to HEK, during 
this stage of the project, given that SC data strongly en- 
hances our capability to detect exomoons. As we men- 
tioned before though, more quarters are becoming public 
over time, including potentially much more SC data, and 
thus the effectiveness of our program is likely to improve 
significantly over time. 

To look for exomoons, we require systems where at 
least three transits have been observed. This is because 
three transits represents the minimum number for timing 
deviations to be inferred. Two transits only introduces a 
total degeneracy between the orbital period of the planet 
and the transit timing variation (TTV) amplitude. Let 
us denote the total baseline of data as B days and assume 
optimistically that Q0 exists. In order to "guarantee'0 
three transits have been recorded, we set the maximum 
allowed orbital period to be Pp < B/4, which is our first 
filter. 

4.4.2. Reliability 

Kepler measures the ratio-of-radii of planetary candi- 
dates, rather than their absolu te radii. From photo met- 
ric measurements discussed in IBrown et al.l (|2011| ) , the 
Kepler input catalogue (KIC) contains estimated masses 
and radii for all of the host stars. Therefore, the plan- 
etary radii can be determined to a reasonable degree of 
reliability. 

However, one of the key parameters in both predicting 
the feasibility of a moon around a planet and the de- 
tectability is the mass of the host planet, Mp (particu- 
larly for computing the Hill sphere). Kepler cannot mea- 
sure the masses of these objects, e xcept in some extre me 
cases where ellipsoidal variat ions ([Mislis et al.l 1201 IT) or 
TTVs exist (lAgol et al.ll2005h . 

The maximum allowed exomoon mass around a planet 

scales as Ms jmax ~ Alp 3 (see Equation [T]) . If we over- 
estimate this value, we will waste time looking at sys- 
tems where moons cannot exist. If we underestimate this 

11 This is not strictly a guarantee due to the potential for data 
gaps but nevertheless serves as a useful minimum requirement. 
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value, we will miss some potentially interesting candi- 
dates. However, the number of candidates is large and we 
can afford to miss some objects for the sake of focussing 
on the very best candidates. Therefore, we choose to 
provide a minimum estimate of Mp. 

Planets are frequently broken up into three regimes: 
i) Super-Earths ii) Neptunes hi) Jupiters. For a Super- 
Earth, the mass-scaling relationship is arguably the most 
predictable relative to the other two regimes, due to the 
known properties of rock under pressure. The precise 
location of the boundary between a gas/ice giant and 
a rocky Super-Earth, and the conditions under which 
such a boundary varies, is empirically a subject in its 
infancy. However, models of the growth of a proto- 
planet's envelope s uggest a critical mass of Mp ~ 20 M® 
(see review article iD'Angelo et al.l[20To| ). whereafter the 
envelope growth timescale rapidly diminishes by orders 
of magnitude. This suggests that the maximum ra- 
dius of a rock y planet is Rp ~ 2i? m , using the scal- 
ing relation of iValencia et al.l (|2006() . This boundary 
is consistent wit h that adopted in other works, such 
as iBorucki et al.l (| 2 1 If ) . Therefore, if a planet has 
Rp < 2i? ffi , we consider it more likely to be icy /rocky 
rather than gaseous. The ass umed mass of t he Su per- 
Earth is calculated using the IValencia et al.l (|2006ft re- 
lation (Rp/R®) ~ (Mp/M e ) - 27 . We stress that this 
assumption is only used for target selection and is not 
adopted in the actual fits or system analyses. 

The definition of a Neptune is somewhat more ten- 
uous but we use the same definition as that of 
IBorucki et al.l ()2011ft - planets which satisfy 2i?® < 
Rp < 6R®. Transiting "Neptunes" are in short sup- 
ply in the exoplanet literature, especially those with 
well-determined densities (SNR> 3)0- Figure [3] shows 
bulk density versus the radius of the six exoplanets 
which sati sfy this criteria. Systems parameters ar e 
taken fromlGillon et all (12011ft. iKundurthv et all (12011ft 



55-Cnc-e 



iFossey et all (12012ft. iBakos et all (12009ft iHe nrv et al.1 
(12011ft and lCochran et al.l (|2011ft for 55-Cnc-e, GJ 1214b, 
GJ 436b, HAT-P-llb, HD 97658b and Kepler-18c respec- 
tively. 

From these, four out of the six settle around an ap- 
proximately common density close to that of Neptune 
and Uranus of pp = 1.7 ± 0.3 gem -3 (where the uncer- 
tainty is the standard deviation of these four). The me- 
dian of all six is 1.5 gem -3 . 55-Cnc-e has a much higher 
density of 4.0^3 g cm ~ 3 , possibly due to its close prox- 
imity (in radius) to the Super-Earth/Neptune boundary 
(Rp — 2.17 ± 0.10 i?®) and is ostensibly not a typical 
member of the 2 — > 6 i?® category. Similarly, Kepler- 
18c is the largest radii planet (Rp = 5.5 ± 0.3,i?®) close 
to the Neptune/ Jupiter boundary of 6i?® and exhibits 
a much lower density of 0.6 ± 0.1 gem -3 , and thus may 
not be typical either. 

We therefore decide to settle on an assumed density 
of 1.7 gem -3 to get an assumed planetary mass in TSA. 
The function of this assumed mass will be primarily in 
computing the maximum allowed exomoon mass and a 
generous margin will be assigned to this calculation to 
account for the fact we have made such a broad approxi- 

12 We define SNR as the reported value for the density divided 
by its uncertainty 
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Fig. 3. — Bulk densities of all known transiting exoplanets in 
the "Neptune" regime. We only show planets with densities deter- 
mined with SNR> 3. Triangles mark the position of Neptune and 
Uranus for context. Vertical grid lines mark th e Supe r-Earth and 
Jupiter boundaries, as defined by Bo rucki et all (201 J ). We mark 
the location of the somewhat anomalous 55-Cnc-e and Kepler-18c. 
The gray ellipse represents a 2-o bound of the bulk of points yield- 
ing an approximately uniform density of 1.7 ± 0.3 g cm~ s . 



mation (as will be discussed later). As before, this again 
only used for TSA and not employed in the final analyses. 

Planets above the 6 i?® boundary show much 
greater variation in density, going from extrema of 
0.122l^gcm- 3 for WASP-17b (jAnderson eU^MHh 
to 26.4 ± 5.6 gem -3 for CoRoT-3b (jDeleuil et al.ll2008ft . 
At this point, we consider the reliability of any mass- 
radius relation to be untenable and thus we add the fil- 
ter than any planet of radius Rp > 6 i?® are not con- 
sidered. This also helps in the exomoon detection too, 
since smaller planets p resent larger m utual events with 
a putative companion (Kipping l2011aft . 

It is important to stress that in all cases that the as- 
sumed masses are only used for target selection purposes 
and will not be adopted as final values in the actual sys- 
tem analyses. 

4.4.3. Capability 

For each Kepler candidate, we need to know the 
capability of the planet to host a moon. The 
maximum allowed exom oon mass from Equation [TJ 
([Barnes fc O'Brienl 120021 ) offers a useful way of accom- 
plishing this. We assume a retrograde moon so that the 
maximum allowed semi-major axis of the moon's orbit is 
93.09% of the Hill radius (i.e. S max = 0.9309). 

To compute the maximum allowed m oon mass, we 
use th e system parameters presented in IBorucki et al.l 
(|2011ft or more recent results where available. As dis- 
cussed earlier, the planetary mass is assumed following 
some simple scaling rules. Finally, we assume the same 
Jovian-like tidal d is sipati on parameters as that used by 
iBarnes fc O'Brienl (|2002ft : specifically Q P = 10 5 and 
k2n = 0-51- Th e fop value is that for a n = 1 polytrope 
(Hubbard 1984) an d the Q p factor is consist ent with es- 
timates for Jupiter (jGoldreich fc Soterlll966ft . This then 
allows us to compute the maximum allowed exomoon 
mass, assuming TT = 5 Gy r, using Equation [TJ taken from 
IBarnes fc O 'Brienl ((20021) . 

It is important to realize this calculation is intended to 
only point the way towards potentially interesting tar- 
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gets. It is not intended to be the final determination of 
this value, which is simply not possible with the current 
information available. 

Accordingly, we exclude all candidates for which the 
maximum moon mass is below 10 M®. This high limit 
is deliberately an overestimate to allow for the number 
of assumptions and approximations made thus far. We 
further exclude planets with clb* < 0.1 AU base d upon 
the arguments given in ^2.5l and ()Namounill2010D . 

4.4.4. Detectability 

In order to make an exomoon detection, we require 
a mass and radius. The timing amplitudes, which yield 
the exomoon mass, depend upon the configuration of the 
system but the eclipse amplitude, yielding the exomoon 
radius, has a far weaker correlation to these terms. The 
eclipse amplitude is approximately given by (Rs/R*) 2 
and thus for any given target we should be able to easily 
estimate the eclipse signal of an Earth-sized moon i.e. 

(R®/R*) 2 - 

We may also estimate the SNR over a 6.5-hour integra- 
tion time by simply dividing the e xomoon eclips e dept h 
by the CDPP values presented in iBorucki et al. (2011). 
We filter all out results where SNR< 1 for a single event. 

4.5. Opportunity Target Selection (TSO) Overview 

Additionally, we will consider "targets of opportunity" 
for special objects of interest. We envision these will typ- 
ically be confirmed, published Kepler exoplanets. These 
targets offer numerous advantages in that the entire pho- 
tometric time series used in the discovery paper is usu- 
ally available (i.e. many more quarters than normal), SC 
data is often available, the planet is known to be gen- 
uine and frequently follow-up information such as spec- 
troscopy, radial velocities or even asteroseismology are 
usually available too. We envision that these TSO tar- 
gets would be typically selected for HEK analysis because 
they have special significance (e.g. are in the habitable- 
zone). 

4.6. Target Selection Prioritization (TSP) Overview 

The prioritization stage is the process of selecting just a 
few targets out of the candidates found by the TSA, TSV 
and TSO stages. TSP is typically done using detailed 
light curve inspection, experience of what constitutes a 
viable signal and other factors. 

For example the availability of short-cadence (SC) data 
is a key TSP factor, due to the improved sensitivity rela- 
tive to long-cadence (LC) data. This is because many of 
the exomoon features induced on the light curve can oc- 
cur on timescales shorter than 30 minutes and thus would 
be lost in the LC data. Further, SC data yields higher 
resolution of the ingress/egress of the planetary transit 
which thus y ields a tighter d etermination of the planetary 
parameters ([Ri pping 2010b). With lower uncertainty on 
the planetary signal, we are naturally able to more eas- 
ily distinguish exomoon signals. Therefore, targets with 
even partial SC data are strongly preferred to those with 
exclusively LC photometry. 

5. FITTING 
5.1. Modelling Strategy: LUNA 



Modelling the eclipses of planet-moon systems is non- 
trivial, if one wishes to retain analytic expressions. The 
advantages of an analytic model are manifold, allow- 
ing CPU intensive fitting techniques (e.g. Monte Carlo 
methods) to fully explore the complex parameter space. 
The requirement for an analytic a l gorithm excludes the 
methods presented in iSimon et al.l (120 09), Sat ofc Asadal 
(2009) , iDeegl (|2009h , and lTusnski fc Valid (12011 7" 

The most significant challenge, in terms of analytic 
modelling, is when the planet, moon and star all partially 
overlap. The analytic solution for the are a of overlap o f 
three circles was only recently found by iFewelll (2006) 
and this discovery allowed for the first time an exomoon 
c ode which could be totally analytic in nature. 

Kipping (|2011aD presented an algorithm to this end, 
dubbed LUNA, which dyna mically mode ls the planet- 
moon motion and utilizes the lFewell (l2006fl solution (p lus 
numerous new solutions derived in iKippina 1201 laf ) to 
produce simulated light curves for moons. The code 
uses quadratic limb darkening and runs almost as fast 
as generating a plane t signal by itself (i.e. the code of 
iMandel fc Agoll 120021 ) . As a result, LUNA is easily im- 
plemented with Monte Carlo based fitting techniques. 
Further, the dynamical component of LUNA means that 
effects such as TTV, TDV-TIP and TDV-V are all inher- 
ently accounted for, plus other previously unconsidered 
effects such as ingress/egress asymmetry. LUNA is a po- 
tent weapon in exomoon detection. 

We note that another analytic algorithm capabl e of 
model ling exomoon eclipses has appeared recently in lPall 
(|2011l) . However, the HEK project will make use of LUNA 
alone, since this already satisfies all of our requirements. 

LUNA also features several other light curve analysis 
techniques developed recently, such as accounting for 
the finite integ ration time using selective resampling 
(Kipping 2010b) and acc ounting for blended / third - light 
using the methodology of iKipping fc Tinettil (|2010D . 

5.2. Detection Strategy: Bayesian Model Selection 

The process of making a detection of any physical phe- 
nomenon is essentially an exercise in model selection. In 
our case, this is most simply described by comparing how 
well the data are explained by a planet-only model (the 
null hypothesis) versus a planet-with-moon model. 

The Bayesian framework is a very powerful basis for 
these model comparisons, allowing the observer to in- 
corporate prior knowledge (such as the allowed phys- 
ical bounds of various parameters) and naturally in- 
clude "Occam's razor" as a way of penalizing overly- 
complicated models. Whilst various information crite- 
rion have been proposed for performing model selection, 
the use of Bayesian evidence has emerged as the met- 
ric of choice to perform model comparisons (jLiddle et al.l 
l200l . 

Bayesian inference methods provide a consistent ap- 
proach to the estimation of a set parameters in a model 
M. for the data D. Bayes' theorem states that 

■Me|D,*,. g-^glM . (14) 

where Pr(0|D, M) — V{&) is the posterior probabil- 
ity distribution of the parameters, Pr(D|0, M) — £(©) 
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is the likelihood, Pr(0|A / J) = tt(0) is the prior, and 
Pr(D|.M) = Z is the Bayesian evidence. 

The evidence can be understood to be simply the factor 
required to normalize the posterior over 0, so that 

Z = J £(0)7r(0)d D 0, (15) 

where D is the dimensionality of the parameter space. 
Since the evidence is the average of the likelihood over the 
prior, Occam's razor is inherently included. Therefore, 
a simpler theory with more compact parameter space 
will yield a larger evidence than a more intricate the- 
ory. Model selection between two competing theories, 
Mo and M.\ can be decided by comparing their respec- 
tive posterior probabilities, for the given data, via 

PrCMijD) gi PrjMj 

Pr(Mo|D) Z Pr(.Mo)' [ ' 

where Pr(A^i)/Pr(A^o) is the a-priori probability ratio 
for the two models, typically set to unity but occasionally 
requires more thought. In this way, the odds ratio of a 
planet-with-moon model can be assessed, relative to a 
planet-only model. Defining the Bayes' factor as B = 
Zi/Zq, I logS| > 6 indicates a > 3 a detection, | log SI > 
10 indicates > 4tr and | logZ?| > 15 indicates > 5cr. 

5.3. Fitting Strategy: MultiNest 

In the exoplanet literature, Markov Chain Monte Carlo 
(MCMC) techniques have emerged as the favoured tool 
for fitting transit light curves. However, the tech- 
nique is primarily for parameter estimation and does 
not natively yield the Bayesian evidence. Whilst tech- 
niques such as thermodynamic integration (e.g. see 
10 Ruanaidh fc Fitzgeraldll201lD can get round this prob- 
lem, it comes at the cost of great computational expense. 

5.3.1. Nested sampling 

Nested sampling (Skillingj I2004J ) is a Monte Carlo 
method which puts the calculation of the Bayesian ev- 
idence in a central role, but also produces posterior in- 
ferences as a by-product. Nested sampling is generally 
considerably more efficient than MCM C methods. For 
examp le, in cosmological applications, iMukheriee et al.l 
(2006) showed that their implementation of the method 
requires a factor of ~ 100 fewer posterior evaluations 
than thermodynamic integration. 

A fu ll disc ussion of nes t ed sam pling is given in lSkillingl 
(|2004t l and IFeroz et al.l il200!);i :. We her e provide a 
brief d escription, following the notation of IFeroz et al.l 
(|2009aft . for the purposes of conceptually illustrating the 
technique. 

Nested sampling takes advantage of the relationship 
between the likelihood and the prior volume to transform 
the multidimensional evidence integral (Equation I15[) 
into a more manageable one-dimensional integral. The 
"prior volume" X is defined by dX = n(<d)d D <d, such 
that 

X(X) = f tt(0) d D ®, (17) 

JC(®)>\ 




(a) (b) 



Fig. 4. — Cartoon illustrating (a) the posterior of a two dimen- 
sional problem; and (b) the transformed C(X) function where the 
prior volumes Xi are associated with each likelihood 

where the integral extends over the region(s) of pa- 
rameter space contained within the iso- likelihood contour 
£(0) = A. One may then re-write the evidence integral 
of Equation [15] in the form 

Z = [ £(X) dX, (18) 
Jo 

where C(X), which is the inverse of Equation [17\ is 
a monotonically decreasing and continuous function of 
X. Consequently, if one can evaluate a series of M like- 
lihoods via d — C(Xi), where Xi is a sequence of de- 
creasing volumes from unity (Xq) to zero (Xm) as shown 
schematically in Figure [4j then the evidence may be ap- 
proximated numerically as a weighted sum using 

M 

Z = y £ j C i w i . (19) 

i=l 

The weights may be computed using the trapezium 
rule, Wi = ^(Xi-i — Xi + \). The algorithm works by 
casting a net of N "active" points across the initial prior 
space. The active point with the lowest likelihood is re- 
moved (made "inactive") and a new replacement point 
is generated such that its likelihood is higher than this 
rejected value. As the algorithm progresses, it travels 
through nested shells of likelihood as the prior volume 
is reduced. The routine is terminated once the evidence 
is computed to some tolerance precision, typically 0.5 
in log-evidence. Upon termination, parameter posteri- 
ors may also be computed using the active and inactive 
points. 

5.3.2. Multimodal nested sampling 

Multimodal nested sampling is an implementation of 
nested sampling to both account for multiple modes 
and achieve efficient sampling via the use of the "si- 
multaneous ellipso idal sampling" method, described in 
IFeroz et al.l (|2007t ). A publicly available version of the 
algorithm, d ubbed MultiNest, is available. We direct 
the reader to IFeroz et al.l (|2009afi for a description of the 
algorithm and its use (including several toy examples) . 

To date, applications of the technique have been mostly 
limited to cosmology, gr avitational wave detection and 
particl e ph ysics (e.g. seelVegetti et al.ll2010l IFeroz et all 
I2009bl and lAbdussalam et all l2(H(T respectivelv) . Re- 
cently however . IFeroz et all (|2011|) demonstrated the first 
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application of the technique to radial velocity data for de- 
tecting extrasolar planets. Here, we briefly discuss our 
application of MultiNest to transit light curve fitting 
(we note that a detailed study of MultiNest for fit- 
ting transits is currently in preparation by Balan et al., 
personal communication). 

5.3.3. Combining MultiNest and LUNA 

Nested sampling allows one to compute the Bayesian 
evidence of a model fit at a much lower computational 
cost than using thermodynamic integration with MCMC 
techniques. However, there are many other advantages 
of employing MultiNest with LUNA, rather than an 
MCMC. 

The most simple and common flavour of MCMC is 
the Metropolis-Hastings algorithm, widely used in the 
recent exoplanet literature. With a unimodal likelihood 
function, this technique is effective and robust, capable 
of identifying the single minimum even when the initial 
starting point of the chain is widely separated from this 
minimum. However, multimodal distributions are ex- 
tremely problematic. If the spacing between two modes 
is much greater than the width of the proposal distribu- 
tion, then the MCMC will take an inordinate amount of 
time to cross over and practically speaking the MCMC 
is stuck in a local minimum. For a planet-only model, 
a unimodal likelihood distribution is generally expected 
but a planet-with-moon model exhibits many modes, es- 
pecially due to harmonic power in the TTVs and TDVs 
(|KiDDind[2009al) . 

MultiNest comprehensively searches the entire prior 
volume identifying all minima and thus is not affected 
by this problem. We point out, however, that more elab- 
orate flavours of MCMC can still find the global mini- 
mum but at further co mputational cost. Techniques such 
as parallel tempe ring ([Gcvcr 1991), differ ential evolution 
(Ter Braak 200(f) and genetic crossovers (|Gregory||2009f) 
can be appended into the MCMC methodology. 

To implement LUNA with MultiNest, we simply need 
to define a likelihood function. MultiNest calls LUNA 
for each active point and LUNA returns a likelihood value 
of the point. It is only at this point where communi- 
cation between the two algorithms occurs and so is the 
only outstanding problem in implementing a combina- 
tion of MultiNest with LUNA. For quiet stars, the noise 
prop erties of the photometry is nearly perfectly Gaus- 
sian ([Kipping Sz Bakos l2011bf l. For Gaussian uncorre- 
cted noise, the likelihood function may be simply defined 
as 

rrr^\ I I ^ \ (/obs,i /mod,i(©)) 

^"'-n^-i a? — J> 

(20) 

where / bs,i is the observed flux, / mo d,i(®) is the model 
flux returned by LUNA, and Oi is the photometric uncer- 
tainty. 

In cases where correlated noise constitutes a signifi- 
cant component of the noise budget, then the likelihood 
function may be m odified accordingly. For example, 
ICarter et al.l (|2008T ) present a wavelet-based likelihood 
function to compensate for time-correlated noise. 



5.4. Choosing the Parameters 

The parameter set for the fitting procedure should be 
physically-motivated so that well-defined boundary con- 
ditions may be imposed, have a uniform prior and exhibit 
as small as possible mutual correlations with the other 
fitting terms. Unless the algorithm is imposed with non- 
uniform priors, the choice of parameters often defines the 
choice of priors too. 

5.4.1. Planet parameters 

The choice of the parameter set for fitting a planet- 
only model is a subject which has been investigated 
in numerous papers in t he exoplanet lit erature (e.g. 
see ICarter et all 120081 and IKippinsj I2010aj) . The tran- 
sit of a spherical planet on a circular orbit across a 
uniformly bright, spherical star is described by just 
four parameters, which is expected given the near- 
trapezoidal nature of the light curve. These parame- 
ters are {tb*,p, (o,b*/R*), where tb* is the time 
of the transit minimum of the planet's barycentre across 
the star, p is the ratio-of-radii (i.e. i?p/i?*), (as*/R*) 
is the orbital semi- major axis of the planet's barycen- 
tre around the star in units of the stellar radius and 
t>B* = (as*/ R*) cosib* (where %b* is the orbital inclina- 
tion of the planet's barycentre around the star). Mul- 
tiple transits allows one to fit for the orbital period, 
Pb*, as well. Of these terms, only (as*/R*) and bs* 
are problematic, due to their very strong mutual correla- 
tion. Further, whilst bs* is bound to be < bs* < 1 for 
fully transiting planets, (aB*/R*) lies within the range 
< (a,B*/R*) < oo and is therefore not appropriate for 
a uniformly dis t ribut ed prior. 

ICarter et al.l (|2008[ ) have suggested using transit du- 
ration related replacements for (as*/R*) and bs* to 
reduce the mutual correlation, b ut these also h a ve es - 
sentially unbounded upper limits. iKipping et al.l (|2012f ) 
have suggested replacing (ob*/R*) with [p" rc ] 2 / 3 (the 
mean stellar density of a spherical star assuming a cir- 
cular orbit, to the power of two-thirds). Besides from 
reducing the correlation to bs*, the bounded limits are 
far better known if we assume the planet is orbiting a 
main-sequence star and impose an upper limit on the ec- 
centricity. Further, the posteriors may be used directly 
to conduct multibody asteroseismology profiling (MAP) 
analysis, which is useful in constrain ing eccentricity in 
multiple systems (Kipping et al. 2012). Thus, we choose 

to use {r B *,p, [p" rc ] 2/3 ,&B*,-Ps*}, and adopt 

irc _ 37r[(aW^) circ ] 3 , 91 . 
— GP B ' 

where we have assumed Mp <C M* to simplify the 
expression. To allow for non-circular orbits, we fit for 
y/es* sin ujb* and ^e B * cosuib*, which maintain uniform 
priors in e B * and lob* but generally exhibit lower mu- 
tual correlations than fitting for e and w directly. Limb 
darkening may be included by adopting a quadratic limb 
darkening law characterized by ui and U2 limb darken- 
ing coefficients. These are constrai ned to lie within th e 
range < u\ + U2 < 1 and u± > (|Carter et alj [2009). 
An upper bound on ui may be estimated by inspection 
of a typical set of coe fficients from stellar atmopshere 
models. For example, iClaretl (|2000D have a maximum 



I. The Hunt for Exomoons with Kepler (HEK) 



13 



u± coefficient of 1.4336 across all listed T e g, logg, etc 
values and across all bandpasses. We use u\ < 2 as a 
conservative upper limit. 

Finally, we allow each transit epoch to have its 
own unique out-of-transit baseline normalization factor, 
OOT. These OOT values are represented by the vector 
OOT are typically not presented in final results tables, 
but are available upon request. 

5.4.2. Moon parameters 

For a planet-with-moon model, a similar set of param- 
eter arises, but with subtle differences. For example, the 
inclination of the moon around the planet- moon barycen- 
tre, isB, can be redefined in terms of an impact param- 
eter as before (e.g. b$B = {asB/R*)cosisB), but the 
< bsB < 1 bound no longer applies, rather we have 
— oo < bsB < oo. Inclination follows < %sb < n 
for prograde moons and u < %sb < 27r for retrograde 
moons. Since neither arcsin nor arccos are unique over 
the range to 2ir and the light curve is not symmetric at 
any boundary, then neither is appropriate for exomoons. 
Instead, we simply use isB- We may also use the full 
range of to 2tt for the bounds to allow both prograde 
and retrograde solutions to be sought. The same is true 
for the longitude of the ascending node, £Isb, which is 
bound by — tt/2 < Hsb < tt/2. Allowing for inclined 
moons is important since an exoplanet's obliquity is not 
expected to tidally decay excep t for orbits relatively close 
to the star (iHeller et al.ll201l( l. 

Eccentricity terms such as esB and ujsb may be fitted 
either directly or using the replacements y'esB sin ujsb 
and iJesB cos ujsb (which still maintains a uniform prior 
in both terms). The latter option tends to be less corre- 
lated and thus preferable. 

Psb may be fitted either directly or fitting for log Psb 
to impose a Jeffrey's prior. In either case, the parameter 
may range from a moon grazing the planetary surface 
(~ hours timescale) to being at exactly one Hill radius, 
occurring at Psb — Pb*/V3 (Kipping 2009a). The phase 
angle of the moon, <fisB, is simply fitted directly within 
the range of to 2n. {cisb/R*) raises similar problems 
as was encountered with facing (clb*/R*)- However, we 
may again make use of the same density trick; namely, 
through Kepler's Third Law we have 



where we have assumed Ms -C Mp to simplify the 
expression. Physically motivated and sensible bounds 
on pp may be easily estimated, in addition to Psb and 
p. It is convenient to fit for the mass and radius ratios 
directly using Ms/Mp and Rs/Rp, since the boundary 
conditions are given by zero to unity in both cases. Since 
the mass and radii ratios are related to the density ratio, 
we choose to keep the planetary density analogous to 
the stellar density by putting to the power of two-thirds 

again i.e. we fit for p 2 J 3 . Table Q] summarizes the priors 
used. 

6. VETTING 

There are two principal signals we are trying to detect: 
timing variations and eclipse effects. Timing variations 



TABLE 1 

Planet-moon parameters used in light curve fits and their 
associated priors. 



Parameter 


Prior 


Planet Parameters 


P 


W{0,0.25} 


[ p circ]2/3 [ kg 2/3 m -2] 


U{ [p* irc ' min ] 2 / 3 , [p» irc,max ] 2 / 3 } 


Ob* 


«{o',i} 


Pb* [days] 




r Bt [BJD TDB ] 


7/JVmin maxT 
U \ T B» ' t B* i 


y/es* coso; Sh , 


«{-l,l} 


yJe B * sin Lois. 




(«1 + U2) 


U{0,1} 


Ml 


«{0,2} 


OOT 


W{0.95, 1.05} 


Moon Parameters 


Rs/Rp 


U{0,1} 


M s IMp 


U{0,1} 


[p P ] 2 / 3 [kg 2 / 3 m- 2 ] 


U ^mi^2/ 3j [p»»]2/3} 


i S B [rads] 


«{0, 2tt} 


U S B [rads] 


W{-7t/2,7t/2} 


Psb [days] 


W{0.083, P B */V3} 


<j> SB [rads] 


U{0, 2tt} 


V e SB cos uj sb 




V e SB sin a; S b 


«{-!,!} 



which could mimic moons can be induced by other per- 
turbing bodies or even stellar activity. Similarly, light 
curve distortions similar in morphology to mutual events 
can be induced by star spot crossings. However, it should 
be noted that out-of-transit companion eclipses are much 
harder to mimic. 

Nevertheless, there is clearly a need to vet exomoon 
signals as being genuine or not, in the same way candi- 
date planets must be vetted. There are numerous tools 
at our disposal to aid in this procedur e. Most g e nerally , 
we follow the six detection criteria of iKippingl (|2011b) 
(see Chapter 1, §4 for details): 

CI Statistically significant 
C2 Systematic errors dealt with 
C3 Physically plausible claim 
C4 Not a suspicious period(s) 
C5 Consistent instrumentation 

C6 Avoid large systematics (e.g. highly active stars) 

Vetting is essentially the act of proposing alternative 
models which can perhaps equally well, or even better, 
explain the observations. Thus we are considering ad- 
ditional models M.2, A^3, M.a, etc and computing their 
evidences relative to that of A4i, the planet-with-moon 
model. However, we note that in some cases, formally 
computing the evidence of an alternative model may not 
be necessary since it may immediately obvious that the 
observations cannot be caused by the alternative model 
in question. 

6.1. Weighing Tools 

The most powerful tool at our disposal for vetting will 
be the derived densities from the light curve fits. As 
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stated in £|3.3l fc B~Tl exomoon systems allow one to mea- 
sure the bulk density of th e star, planet a nd moon all 
through photometry alone (Kipping 2010cJ). If one also 
has some radial velocity data, then the absolute dimen- 
sions of all three bodies can be determined using Kepler's 
Third Law. 

If the moon is not genuine, but caused by some false- 
positive or even just residual noise in the data, then it is 
highly improbable that these values will come out to be 
remotely physical. Thus, by bounding the prior volume 
to only physically plausible values (as discussed in i)5.4[) . 
we forcibly exclude the vast majority of false positives. 

In most cases, we will not have a determination of the 
radial velocity semi- amplitude, K*, at our disposal due 
to the faintness of the Kepler targets (although we may 
have an upper limit). Even without K*, we still directly 
measure p*, pp and ps- Using the some reasonable es- 
timate for M* and allows for the dimensions of all 
three bodies to be determined. 

6.2. Starspot Crossings 

Starspot crossings can mimic exomoon-like mutual 
events. There are several ways in which they differ 
though. First of all, a starspot crossing tend s to be 
V-shaped (e.g. see iSanchis-Oieda fc Winnll2011l) due to 
the fact that the spots are typically greater than or ap- 
proximately equal to the size of the transiting planet. 
In contrast, exomoon signatures should usually be flat- 
topped, but can in some more non-aligned geometries 
be V-shaped. Nevertheless, a flat-topped mutual event 
would be a strong indicator that the event was due to a 
moon rather than starspot crossing. 

Secondly, starspots have chromatic variations whereas 
a mutual event should not. Thus there is the possibil- 
ity of conducting follow-up observations from the ground 
to interrogate this hypothesis. Thirdly, starspots should 
track across the stellar surface with a speed determined 
by the rotational period of the host star, P* >ro t- This 
term is usually determinable from the photometry alone, 
should the star be sufficiently active. 

Since an auxiliary transit outside of the main tran- 
sit event is much more difficult to mimic through stellar 
activity, these events will undoubtedly be more easily de- 
tected. These events tend to be associated with moons 
on large separations such as a wide-captured retrograde 
moon and thus we anticipate that HEK may have a de- 
tection bias in this regard. 

However, starspots arc unlikely to be both well- 
modelled by an exomoon-fit and produce physically plau- 
sible densities for the star, planet and moon. This is 
because starspots tend to produce very strong eclipse 
features but only weak timing artifacts. For exam- 
ple, HAT-P-llb transits a heavily spotted star ex- 
hibiting eclipse features of u p to l-2mmag in depth 
(|Sanchis-Oieda fe Winn! I2011D bu t timing deviations of 
< 30 seconds (|Deming et aL f l20ll . The absence of sig- 
nificant timing variations would cause the planet-with- 
moon fit to favour an implausibly low density moon. 
Thus, the weighing tools will be a most common test 
for the presence of spots. 

6.3. Transit Timing Analysis 

Transit timing variations (TTVs) are induced in both 
resonant (jAgol et al.ll2005t iHolman fc Murravll2005l) and 



non-resonant systems (Nes vornv fc Beaug e 2010). One 
possible model to explain a candidate moon system 
would be a transiting planet influenced by an unseen per- 
turbing planet, A^2- The evidence of this model may be 
computed by allowing each transit epoch to have its own 
unique time of transit minimum, tb*a, and proceeding 
as before with the planet-only fit. 

Aside from producing an evidence value more 
favourable than the planet- with-moon fit, model M.2 
must correspond to a physically plausible scenario. To in- 
vestigat e this, we will use the TTV- inversion code devel- 
oped bv lNesvornv fc Beaugel (|2010l ). ttvim.f , to explore 
the range of plausible perturbers. If no plausible per- 
turbers can fit the transit times, then model M.i will have 
a negligible prior probability i.e. Pr(A^2) — > 0. Such an 
instance would thus favour the planet-with-moon model 
over the perturber model. If the perturber is both feasi- 
ble and of comparable evidence to the moon-model, then 
one may search the photometric time series for transits 
of the other planet directly. 

The same process may be repeated for other sources 
of TTVs too, such as the model of a Trojan perturber 
(jFord fc Holmanl 120071 ). model For such scenar- 

ios, transits of the Tro jan(s) can also be sought (e.g. 
IKipping fc Bakoil2011aH . 

6.4. Radial Velocities 

Whenever possible, the HEK project will seek to obtain 
radial velocities of the target systems to confirm both the 
planetary and lunar nature of the candidates. Since TSA 
targets only sub-Jupiters on moderate-to-long periods, 
the expected RV amplitudes will be typically < lOm/s. 
Detecting this signal would allow us to both confirm the 
candidate and determine the absolute dimensions of the 
system. However, excluding RV amplitudes above a cer- 
tain threshold is also useful in eliminating false positives. 

Aside from dynamically weighing the system, we antic- 
ipate precise radial velocities will be very useful in locat- 
ing the correct orbital period o f the moon. As discussed 
earlier and in Kipping (2 009al ). the timing signals can 
be fitted by a forest of harmonic orbital periods due to 
the undersampled nature of our observations. A very 
strong eclipse signal can remove this degeneracy, as can 
enforcing coplanarity, but if either of these is not possi- 
ble, we are left with a forest of harmonics. Each period 
yields a unique pp and thus unique light curve derived 
M p/M* . For example, a synthetic example shown later 
in §7.21 finds two harmonic modes with derived pp val- 
ues differing by more than a factor of two. If precision 
radial velocities can infer Mp/M* independently, then 
the correct harmonic is empirically determined. Thus, 
systems with indistinguishable harmonic modes will be 
prioritized for radial velocity follow-up. 

7. WORKED EXAMPLES 

7.1. Null Example: HZ Neptune around an M2 star 

As a null examp le, we use a synt hetic planet-only light 
curve model from Kipping (2011a). The synthetic data 
consists of six SC transits of a habitable-zone (HZ) Nep- 
tune around an M2-star, observed with Kepler-class pho- 
tometry (specifically we used Gaussian noise of 250 ppm 
p er minute ) . Deta ils of the model can be found in §4.6 
of lKippind (|2011al) . 
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We fit the data with two models: Mo, a planet- 
only model and A4i, a planet- with-moon model. In 
both cases we assume circular orbits for simplicity. The 
adopted priors are the same as that given in Table [TJ 
We used [pf in ] 2 / 3 = 18.4982kg 2 / 3 m~ 2 , correspond- 
ing to the lower limit on th e lowest density exo planet 
presently known (WASP-17b: IAnderson et al.|[20Toh . and 
[pp ax ] 2 / 3 = 920.9 76 kg 2 / 3 m~ 2 , corresp onding to a iron- 
rich Super-Earth ([Valencia et al.1 120061 ). For the st ellar 
density, we use the range provided in iCoxl (|2000f ) for 
main-sequence stars of spectral type M5 to F0. Pb* and 
tb* are given priors of ±1 d around the known solution. 
The solution for these terms is always easily inferred from 
simple inspection of a photometric time series. 

The Bayesian evidence for the two models was found to 
be log(Zo) = 23377.36 ± 0.21 and log(-Zi) = 23374.76 ± 
0.22, thus giving |A[logS]| = 2.60 ± 0.30 in favour of 
the planet-only model. This occurs despite the fact the 
planet-with-moon model yields a lower % 2 of 8921.57, 
versus 8930.96 for the planet-only model, thus demon- 
strating the built-in Occam's razor of Bayesian evidence. 

Despite the Bayesian evidence clearly favouring the 
planet-only model, as expected, the posteriors from the 
planet-with-moon model may be used to place upper 
limits on the allowed mass and radius of the exomoon. 
The posteriors, shown in Figure [6j exclude M$/Mp > 
0.018 and R s /Rp > 0.44 to 90% confidence (the 4- 
<7 constraints are of limited use; Mg/Mp > 0.369 and 
Rs/Rp > 0.999). 

7.2. Moon Example: HZ Neptune around an M2 star 
with a Widely Separated Earth-like Moon 

As a moon example, we u se a synthetic p lanet-with- 
moon light curve model from iKippind (|2011al ). The syn- 
thetic data consists of six SC transits of a habitable- 
zone (HZ) Neptune around an M2-star (as before), ob- 
served with Kepler-class photometry (same noise as be- 
fore). The moon is set to an Earth- mass/radius object at 
the edge of the Hill sphere on a retr ograde o r bit. De tails 
of the model can be found in §4.5 of lKippina (|2011a[ ) and 
the de-noised model is shown in the left panel of Figure[2l 

Fitting the data using MultiNest with a planet- 
only model yields a Bayesian evidence of logZ = 
22104.55 ± 0.21, whereas the planet-with-moon model 
reaps logZi = 23552.75 ± 0.27. The corresponding 
posteriors of both models are shown in Figures EMHl 
The change in evidence corresponds to a 54-ct detection, 
which is clo se to the 50-g d etection found using a sim- 
ple F-test in IKipping] (|2011af) . However, MultiNest re- 
veals that three distinct modes exist in the data, which 
are reported in Table [2j Modes 1 & 3 correspond the 
correct exomoon orbital period of 23.995 d, but mode 2 
is located at Psb — 15.775 d. Th is corresponds to a har- 
monic, as originally predicted in IKipping (2009a). This 
is confirmed by evaluating the expected position of the 
first harmonic using [{I /Psb) + (l/^s*)] -1 = 15.769d. 

Modes 1 and 3 both locate the true period and both 
exhibit a significantly improved evidence value. In fact, 
mode 2 is disfavoured at the >7-a level over modes 1 and 
3. Between modes 1 and 3, the difference is that mode 3 
occurs at the correct retrograde solution whereas mode 
1 is prograde. The evidence difference between these 
modes is |A(logZ) = 1.2 ± 0.6 i.e. insignificant. In a 



blind analysis, there would be no way of reliably distin- 
guishing the correct solution based upon the available 
data, but more transits should help since the TDV-TIP 
effect is asymmetric with respect to the orbital sense of 
motion. 

Despite the presence of three modes, MultiNest pro- 
vides an appropriately weighted global posterior as well 
as the local posteriors. Using this weighted posterior 
heavily deweights mode 2 and the leads to accurate val- 
ues for the key terms such as mass, radius and density 
(see last column of Table [2] Orbital angles such as incli- 
nation become unconstrained due to mixing the prograde 
and retrograde solutions together, but this is an accurate 
representation of our ignorance of this term. 

8. SUMMARY 

In this work, we have presented a description of the ob- 
jectives and methods of our new observational project, 
"The Hunt for Exomoons with Kepler" (HEK). The HEK 
project will seek to infer the presence of extrasolar moons 
around transiting exoplanet candidates observed by the 
Kepler Mission. In cases of null detections, upper lim- 
its will be reported and the full set of parameter pos- 
teriors will be made available on the project website 
(www.cfa.harvard.edu/HEK/). These statistics may be 
used to deduce the frequency of large moons around vi- 
able exoplanet hosts, r/^ . 

We have described, in §3J how the list of all Kepler 
Objects of Interest (KOIs) will be distilled into a subset 
of the most promising candidates, for the purposes of ex- 
omoon detection, via a three-prong target selection (TS) 
strategy. This includes visual identification (TSV), auto- 
matic filtering (TSA) and targets-of-opportunity (TSO) 
with a final stage of target prioritization (TSP). 

Selected targets will be interrogated for evidence of an 
exomoon by comparing the Bayesian evidence, Z of a 
planet-only model and a planet-with-moon model (see 
discussion in q5.2[) . In addition to presenting a > 4- 
o~ preference for the planet-with-moon model, putative 
candidates must demonstrate physically plausible solu- 
tions. 

In fitting the data, we require both a forward-model 
and a fitting algorithm. The forme r task is handle d 
by the LUNA algorithm, developed by IKipping (2011a), 
which is designed to analytically model the transit light 
curve of a planet-with-moon system including limb dark- 
ening, dynamical motion and mutual events. Due to 
the highly complex parameter space, featuring multiple 
modes due to aliased harmonic power, and the need for 
an highly efficient computation of the Bayesian evidence, 
a sophisticated fitting algorithm is required. To this 
end, we have pres ented the application of MultiNest 
(jFeroz et al.1 120071 ) with LUNA. MultiNest is a multi- 
modal nested sampling algorithm widely used in the cos- 
mology and particle physics communities. £15.31 describes 
our implementation and the choice of priors, parameter- 
ization and likelihood function. 

Strategics for vetting potential candidates are dis- 
cussed in §32 before we present two examples of 
LUNA+MultiNest on synthetic data in §3 The example 
fits demonstrate not only the multimodal nature of look- 
ing for exomoons but also how Bayesian evidence may 
be used to detect such systems. 

We are currently analyzing a subset of preferred can- 
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didates and will be reporting on these findings in the 
near-future (Kipping et al. 2012, in prep.). As the HEK 
project progresses, we hope to answer the question as to 
whether large moons, possibly even Earth-like habitable 
moons, are common in the Galaxy or not. Enabled by 
the exquisite photometry of Kepler, cxomoons may soon 
move from theoretical musings to objects of empricial 
investigation. 
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TABLE 2 

Comparison of exomoon parameter estimates from three modes found in the MULTLNbst fits of synthetic data. Data generated for a 
Neptune with a distant moon around an M2 star. The global evidence of the planet-with-moon model (all three modes) is log = . Mode 
is the most accurate mode but the local evidence values between modes 1 and 3 are sufficiently close that distinguishing these modes 

blindly would not be possible (mode 2, however, can be discounted). 
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Fig. 7. — Marginalised posteriors from MultiNest when fitting a synthetic example of a planet-with-moon data set using a planet-only 
model from LUNA. 
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Fig. 8. — Marginalised posteriors from MultiNest when fitting a synthetic example of a planet-with-moon data set using a planet-with- 
moon model from LUNA. We here only show the exomoon- related parameters for brevity. 



I. The Hunt for Exomoons with Kepler (HEK) 
APPENDIX 



TABLE 3 List of important acromyms used in this paper. 

Acromyn Definition 

~~ FTEK Hunt tor Exomoons with Kepler project 

SC short-cadence 

LC long-cadence 

KIC Kepler Input Catalogue 

TTV Transit timing variations 

tf$-v teryta TW ions 

TDV-TIP Transit impact parameter induced TDV 

TS Target selection 

TSA Automatic target selection 

TSV Visual target selection 

TSP Prioritization target selection 



TABLE 4 List of important parameters used in this paper. 



Parameter 

ST, 

VSB 
es* 
ess 
wb, 

USB 

«S* 

tSB 

it* 

Rp 
Rs 

a 

a-B* 
a-SB 
P B * 

k B 

M P 

M S 

p* 

Pp 

PS 

V 

s 

6 

TB* 

ob* 
f 

<5ttv 

<5tdv-v 

<5tdv-tip 

Attv 

Atdv-v 

Atdv-tip 

$TTV 

'I'tdv-v 

'I'tdv-tip 

Qp 

k2 P 

T 



Definition 



Sky-projcctcd separation between the planet-moon baryecntre and the host star, in units ot K* 
True anomaly of the satellite around the planet-moon barycentre 

Orbital eccentricity of the barycentre of the planet (+ any satellites) around the host star 
Orbital eccentricity of the satellite around the planet-moon baryecntre 

Argument of periapsis of the baryecntre of the planet (+ any satellites) around the host star 

Argument of periapsis of the satellite around the planet-moon barycentre 

Orbital inclination of the barycentre of the planet (+ any satellites) around the host star 

Orbital inclination of the satellite around the planet-moon baryecntre 

Longitude of the ascending node of a satellite, relative to the orbital plane of the planet 

Radius of the host star 

Radius of the planet 

Radius of the satellite 

Semi-major axis of an object around its primary 

Semi-major axis of the planet-moon barycentre around the host star 

Semi-major axis of the satellite around the planet-moon barycentre 

Orbital period of the baryecntre of the planet (+ any satellites) around the host star 

Orbital period of the satellite around the planet-moon barycentre 

Mass of the star 

Mass of the planet 

Mass of the satellite 

Mean density of the star 

Mean density of the planet 

Mean density of the satellite 

Ratio of the planet's radius to the stellar radius (Rp/R t ) 
Ratio of the satellite's radius to the stellar radius (Rs/R*) 
Defined as p 2 

Instant when dSpt/dt = near inferior conjunction 

Impact parameter of the barycentre of the planet (+ any satellites) 

Time between the planet's centre crossing the stellar limb to exiting under the same condition 

RMS amplitude of TTV signal 

RMS amplitude of TDV-V signal 

RMS amplitude of TDV-TIP signal 

Waveform of TTV signal 

Waveform of TDV-V signal 

Waveform of TDV-TIP signal 

Enhancement factor of TTV signal 

Enhancement factor of TDV-V signal 

Enhancement factor of TDV-TIP signal 

Tidal quality factor of the planet 

Love number of the planet 

Lifetime of the moon 



